Structural phase transition at high temperatures 
in sohd molecular hydrogen and deuterium 
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We study the effect of temperature up to lOOOK on the structure of dense molecular para- 
hydrogen (P-H2) and ori/io-deuterium (0-D2), using the path-integral Monte Carlo method. We find 
a structural phase transition from orientationally disordered hexagonal close packed (hep) to an 
orthorhombic structure of Cmca symmetry before melting. The transition is basically induced by 
thermal fluctuations, but quantum fluctuations of protons (deuterons) are important in determining 
the transition temperature through effectively hardening the intermolecular interaction. We estimate 
the phase line between hep and Cmca phases as well as the melting line of the Cmca solid, (accepted 
for publication in Physical Review B) 
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I. INTRODUCTION 



The possibility for hydrogen to undergo a pressure- 
induced phase transition from a proton-paired insulator 
to a monatomic metal was first suggested in 1935.a Since 
then a large body of experimental and theoretical works 
have been done to determine its phase diagram in the 
pressure (P) and temperature (T) plane. |I)espite studies 
at P reported to be as high as 342 GPa,H the tenacious 
covalent-bond feature of a pair of protons does not allow 
observation of a monatomic phase so far. 

For P up to ca. 200GPa, the phase diagram has 
been investigated by several research groups at the room 
temperature and below, based on the optical measure- 
ment performed in the diamond anvil cell (DAC) de- 
vices. By now it is well established that the solid hy- 
drogen exhibits at least three different molecular phases 
(the phagLes I— III), although some details are still in 
dispute.|jl3 (I) Phase I: At P < IIOGPa, each center of 
the II2 molecule occupies the lattice site of an hep struc- 
ture, but quantum rotational effects overcome librational 
barriers, leading to an orientationally disordered phase. 
(2) Phase II or broken-symmetry phase (BSP): At P be- 
tween llOGPa and 150GPa, anisotropic intermolecular 
interactions freeze the molecular rotations into an orien- 
tationally ordered phase. (3) Phase III or H-A phase: At 
P above 150 GPa, a third phase is observed, expected to 
be another kind of orientationally ordered phase. 

In 1996, Weir and coworkers gave one of the exciting 
results. They observed high electrical conductivity in the 
shock compressed H2 and D2 liquids and interpreted it as 
a transition from a semiconducting to a metallic diatomic 
fluid at 140 GPa and 3000K.y This experiment clearly 
demonstrates the importance of temperature effects in 
the search of metallization. For T larger than 5000K, 
the effects are studied theoretically in the molecular, the 
dissociated, and the plasma regime of dense Il2.E2l't3 It is, 
however, still largely unclear how T influences the state 



of solid H2 at T higher than the room temperature. We 
need to know the phase diagram at this temperature re- 
gion to connect the liquid state with the low-temperature 
phases for a comprehensive understanding of dense H2. 
Thus we focus our attention on the range of P up to 
200GPa and 300K< T < f OOOK in this paper. 

In order to study the condensed II2 phases theoreti- 
cally, several methods have been adopted at various levels 
of approximations to the ab initio Hamiltonian represent- 
ing the coupled system of Na protons and Na electrons. 
These include the calculation of electronic energies in 
the local density approximation (LDA) or its refinements 
to the density functional theory for a variety of crystal 
structures with moleculajLorientations fixed to some par- 
ticular configuration,E3 Ej and the implementation of ab 
initio moie&wdar dynamics treating protons as classical 
particles.cSiHl We should note, however, that the strong 
quantum nature of light protons requires a more care- 
ful quantum-mechanical description of their zero-point 
motion. In fact, using eithjei; the first-principle path- 
integral molecular dynamicsrH. ar the quantum Monte 
Carlo (QMC) simulations,EifEjO the calculation treat- 
ing protons as dynamic quantum particles has already 
been performed. Although desirable, this approach de- 
mands the computer resources very much. Thus, even 
in a recent paper ,0 the calculation is done only with as 
small as Na = 64. This size of Na is too small in order to 
obtain a correct equilibrium structure free from any bias 
of an initially assumed one through simulations. 

A ten times increase in Na is possible if we approach 
the problem by adopting the hydrogen molecule as a 
basic ingredient in simulations rather than the proton- 
electron mixture. In this approach, the system is re- 
duced to a quantum-mechanical problem of N (=2Na) 
molecules interacting to each other through an effec- 
tive intermolecular pair potential Vpair- Thus the QMC 
calculatiouris iHiplemented only for the nuclear degrees 
of freedomE2rEZl and the electronic degrees of freedom 
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are implicitly taken into account in the choice of Vpaii - 
In the present paper, we adopt this approach and em- 
ploy the finite-temperature path-integral Monte Carlo 
(PIMC) method to investigate both lattice and oricnta- 
tional transitions in the molecular phase with zero-point 
motions incorporated rigorously. We shall use an empir- 
ically determined Vpair which is the saipi of the Hemley- 
corrected Silvera-Goldnuui potentiaS and the Runge- 
scaled Shaefer potential. S This choice of Vpair is known 
to give an equation of state in very good agreement with 
experiment as well as the I/II phase boundaries for both 
H2 and D2 over the whole experimentally investigated 
range of pressures .EZI 

This paper is organized as follows: We shall explain 
our theoretical model in more detail in Sec. including 
a review of trials to determine Impair, a description of our 
choice of Ipair, a brief summary of the PIMC method 
with a constant-pressure ensemble, and some computa- 
tional details. In Sec. Ill our results are shown for both 



solid II2 and D2, indicating a temperature-induced solid- 
solid structural phase transition. The structure of the 
new high-temperature phase, the solid-liquid phase tran- 
sition, and the phase diagram are also given here. Section 
[V contains a discussion on the role of zero-point motion 
of protons in the solid-solid structural phase transition 
and finally in Sec. we summarize our results. 



II. THEORETICAL MODEL 

A. Molecule base 

Within the temperature and pressure range that a 
molecule can be regarded as a basic ingredient, a quan- 
tum molecular solid can be described by the Hamiltonian 
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Here is the center-of-mass position vector of the ith 
molecule, fli its orientation vector, its angular mo- 
mentum operator, and the intermolecular separation vec- 
tor R.y = Rj — Rj . The molecular mass and moment of 
inertia are denoted by m and /, respectively. For a lin- 
ear molecule with the intramolecular bond length fixed, 
/ is a scalar constant. The values of m = 3676 and 7352 
(atomic units) for H2 and D2 , respectively, and rotational 
constants B = 121 = 8f£lS,and 42.92 K for H2 and 
D2, respectively, are used.E3Ej Among the the nuclear 
degrees of freedom, we include both translation and ro- 
tation modes of molecules in the kinetic energy Tnucicar, 



but we consider implicitly the intramolecular vibration 
mode along with the electronic degrees of freedom by ju- 
diciously choosing Vpair the pair- wise sum of which gives 
the potential energy Kucioar- 

We have enough evidence to support the present ap- 
proach in our research range (0 ^ 200GPa and ^ 
lOOOK). First, we can neglect the molecular dissociation, 
because its fraction is less than 0.01% from the same 
analysis-as-the one leading to about 5% at 140GPa and 
3000K.t3c3 Secondly, ab initio quantum-chemistry calcu- 
lations at N=2 lead us to the conclusion that the depen- 
dence of V^air on the bond length tq is weak, even though 
To shrinks by more than 10% fcpm 1.4 bohr (0.74A) as 
two II2 molecules come closer .EJ This assures the valid- 
ity to employ an rg-independent Vpair- The notion of 
Vpair is relevant if r the intermolecular distance is larger 
than about 2.3 bohr (I.2A), because each molecule is esti^, 
mated to possess a repulsive core of the radius of 0.6A.EJ 
Finally the effect of zero-point intramolecular vibrations 
(vibrons) can be taken into acc ount in the form of V^pair 
as we shall illustrate in Sec. IV by using a toy model. 



B. H2-H2 interaction 



Rather than determining ypair from first principles. 



we shall resort to a phenomenc 
well reviewed in the literature.! 



deal approach which is 
^1 In general, the multi- 



dimensional function Vpair (R12, S^i, ^^2) can be expressed 
as the sum of the isotropic pair potential Vo(i?i2) and 
the anisotropic potential Va„i(Ri2, ^2). A variety of 
analytic forms have been proposed for Vb(i?) and their 
relevance has been tested in the past. Four of them are il- 
lustrated in Fig. |l|. Because the Lennard-Joncs (LJ) two- 
parameter form (the dotted-dashed curve in Fig. |l|) has 
an unphysically strong repulsive core, more sophisticated 
forms have been proposed. The most, commonly used 
foriTLis either BUCK potential Vbuckc3 or SG potential 
VsG, exploited by dashed and dotted curves, respectively, 
in Fig. ^. Although they have almost the same analytic 
form, Vbuck is fitted to the gas-phase data, while VsG to 
the solid-phase data. In order to include the three-body 
interactions effectively in the solid environment, a fur- 
ther refinement is introduced to the original form of Vsg 
by incorporating an additional repulsive term, leading to 
the final form of SG potential. Compared to Vbuck, the 
net effects contributed from the three-body interaction 
are a hardening of Vsg at small R and a slight raising 
of the well depth. Their concrete forms can be found in 
Refs. |3,|9[ and ||. 

This Vsg (R) works well for the solid II2 and D2 under 
ambient pressures. Under high pressures, however, an ad- 
ditional softening effect found by Hemley et al. is needed 
to correctly describe the enhanced many-body contribu- 
tions in the short-range region caused by the dense solid 
environment. This effect can be incorporatecL-by an ad 
hoc short-range correction VsR(i?) to Vsg(^),^ leading 
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to the definition of the Hemley-corrected SG potential 
Vsgh(^) (the solid curve in Fig. |l|), where the actual 
form of Vsk{R) can be found in Refs. ^ and |3^. This 
softening correction at small R is very important in re- 
producing the equation of state (EOS) derived—experi- 
mentally from the liquid D2 shock-wave data,OEJ as as- 
sessed by Cui et al. in Ref. ^ in which the EOS for 
D2 at T = 20K and for P up to 40GPa is calculated for 
testing the intermolecular potential. From the test, the 
EOS is proved to be insensitive to the details of T4„i, but 
it is very sensitive to the choice of Vq{R) which controls 
the translational motion of molecules. In fact, Vsgh(-R) 
gives reS|tJ.ts in excellent agreement with the experimen- 
tal dataj2a while Vsg(-R) does not. The same conclusion 
is also drawn at T = 300K from a similar test; for P up 
to 200GPa, Vs,GYi{R) provides the resuhs for the EOSLm 
line with the experimentally obtained fitting formula.EZI 




FIG. 1. Various potentials proposed for the interaction be- 
tween two H2 molecules. The solid, dotted, dashed, and 
dotted-dashed curves represent, respectively, the SGH, SG, 
BUCK, and LJ potentials. 

Only a few works have been done for the assessment 
of Vani- Using the ah initio quantum chemical method 
to evaluate the contributions such as the long-range elec- 
tronic quadrupole-quadrupole interaction and the atom- 
diatom scattering_tj3p3ji, Schaefer et al. have given a 
potential for Vani-c3t£l Comparing the results obtained 
from the electronic-structure calculation in the LDA, 
Runge et al. have found that this Schaefer potential 
Vschaefer IS too rcpulsivc in the dense solid. Thus, in 
order to soften it, they have introduced a scaling factor 
to Vschaefer to provide Vani as Vschacfcr-E^ This fac- 
tor a„ is linear in the nearest neighbor spacing i?„„ to 
correctly describe the dense-solid environment as a„ — 
0.61-K).31(i?nn/^°„-0-5) with i?o„ = 3.789A, the H2 equi- 



librium zero-pressure nearest-neighbor spacing. Both the 
bare Vschacfcr and the scaled one for Vani have been tested 
in Ref. 37 to find that this difference is crucial in re- 
producing the I/II phase line for the solid D2. The 
best agreement with the experimental data for P up 
to 150GPa is achieved by using the scaled anisotropic 
potential combined with Vsgh(-R) for Vo(i?). The bare 
Vschacfcr predicts the transition to occure at much lower 
pressures. This diffetjence was also observed in the fixed- 
lattice PIMC study.y 

For those reasons mentioned above, we shall choose 
the sum of I^sgh {R) and the scaled Schaefer potential as 
We believe that this Vpair is most 



Vpair in this paper. 



reliable for the study of solid molecular H2 and D2 in our 
research range, in particular, for P < 150GPa. 



C. Path Integral Monte Carlo method with constant 
pressure ensemble 



With the Hamiltonian thus provided, we model the 
bulk solid by a simulation cell, which is periodically du- 
plicated in all three spatial dimensions to minimize sur- 
face and finite-size effects. Initial size and geometry of 
the simulation cell are chosen to accommodate a partic- 
ular number of molecules (TV) and hep structure. Our 
calculations are performed mostly on TV =288 and from 
an initial simulation cell determined by two basis vec- 
tors (ap and bp) forming a 60" angle and by the third 
one (cp) perpendicular to them with the length ratio of 
ap : 6p : Cp = 1 : 1 : 4-^/6/9. The packing pattern is 
ABAB... to form the hep lattice structure. Extensive 
testing on iV = 96, 150 and 392 has also been made. 

In studying the effect of heating solid H2 isobarically at 
high pressures, we perform the path-integral Monte Carlo 
calculations with a constant-pressure {NPT) ensemble, 
instead of a simpler constant-volume (NVT) one, to 
avoid a bias of the restrictive cell geometry with a prede- 
termined crystal structure. Implementation of the NPT 
ensemble is achieved by extra independent Metropolis 
moves of three basis vectors ap, bp and Cp, which gen- 
erates a Markov chain of states having a limiting distri- 
bution proportional to exp{—H/T + NlnV) at T with 
the enthalpy H{= PV + Eg), where Eg is the energy ex- 
pectation value (iJnucicar) of the configuration s with s 
representing a set of scaled coordinates. This enables us 
to monitor volume changes and therefore to observe a 
possible first-order phase transition directly. 

In order to avoid the "minus-sign problem" in QMC 
studies inherent in fermions, we confine ourselves to 
studying P-II2 and 0-D2. After 4000 Monte Carlo steps 
for equilibration, statistical averages are collected from 
every second step, to a total of about 4000 data points. 
Because we focus on higher temperatures than the room 
temperature, relatively smaller partition number in the 
imaginary-time axis (M) is enough. At a fixed P, we 
start with simulations at the room temperature at which 
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M is increased up to 16 to ensure the convergence. We 
come to know that M = 8(5) is enough for P = 30GPa 
(higher P) . Once its large enough value is determined, M 
is held fixed during heating at the given P. This allows 
us to use the equilibrium configuration at low tempera- 
ture as a start-up one for higher temperature, resulting 
in continuous heating of our sample. Because other com- 
putation techniques are documented in Ref. |3^, we omit 
further details of them here. 



III. RESULTS 



We observe the equilibrium structure in real space di- 
rectly and also monitor the pair distribution function 
g{R) that represents the conditional probability of find- 
ing other molecules at the distance R from a specified 
molecule at the origin: 

9iR) = ^pi^m,-R)^. (3) 

with p the density of system. This function shows a set 
of well-defined peaks characteristic to the configuration 
of neighboring molecules around the one at the origin. 
The function 0{R) defined in Ref. ^ is used to measure 
the correlation in molecular orientation. 



FIG. 2. Snapshots of the average equilibrium distribution 
of molecules at lOOGPa in real space projected on both (i) xy 
and (ii) xz planes. Row (a): p-H2 at 450K; row (b): p-H2 at 
650K. The structure of each unit cell is shown in (iii) in which 
the size of fluctuations is represented by the painted region. 

The equilibrium distribution of 288 (= 6 x 6 x 8) 
molecules in real space projected on both xy and xz 
planes are shown in Fig. || for P-H2 at lOOGPa and at 
two different temperatures. Although it fiuctuates with 
the standard deviation indicated by the painted region 
in Fig. ||(iii), on the average, each molecule occupies the 
lattice site of hep at 450K, while that of an orthorhombic 
Cmca structure at 650K. 

For detecting this transition quantitatively, we plot 
g{R) at P = lOOGPa with the increase of T in Fig. |. 
Below 560K the second peak in g{R) at i? ~ 2.72A, a 
characteristic feature of hep, is clearly seen. But the 
peak disappears at higher T, which is a sign of the struc- 
tural transition from hep to Cmca. Thus we identify 
the transition temperature Ttr as 560K with a statistical 
error of 20K. (By "statistical" we mean that the differ- 
ence of states cannot be seen clearly if that in T is less 
than 20K due to statistical fluctuations.) This structural 
change is always seen in the system with N = 96, 150, 
or 392. We find a large N dependence of Ttr, leading to 
an error of 40K or larger incurred at extrapolation of Ttr 
at N 00. This error cannot be reduced further, be- 
cause calculations with a much larger A'' are not feasible 
at present. 
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FIG. 3. Temperature dependence of g{R) for P-H2 at 
lOOGPa. Inset shows the temperature dependence of volume. 
(Statistical errors are smaller than the size of symbols.) 
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With the further increase of T from Ttr, the peaks in 
g{R) become less sharp. EventuaUy at T above 875K 
(which we identify the melting temperature with a 
statistical error of 25K, about the same size of an error 
due to extrapolation in A'^ — > c») g{R) exhibits charac- 
teristics of a liquid phase. This first-order melting tran- 
sition can be much better identified by the discontinuous 
change with T in either the enthalpy H or volume V (the 
inset of Fig. ||). 
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FIG. 4. (a)Phase diagram of P-H2. At low temperatures, 
the dashed curve shows the I/II phase boundary interpolated 
through the two measured data marked by open diamonds 
(Refs.^ and while the dotted-dashed curve indicates the 
experimental II/III phase boundary. Solid diamonds repre- 
sent the orientational order-disorder phase transition points 
obtained by PIMC in Ref. ^ (Size of the symbols represents 
the magnitude of errors.) Solid triangles and open circles for 
the hep/ Cmca transition and the melting, respectively, show 
the new phase boundaries obtained in this study. The dot- 
ted curve is a Simon fit to the melting line of this work, (b) 
Meltin&4iiie at low pressures wiA experimental data for n-H.2 
(pluses)Ell and for P-H2 (stars). EJ 

By evaluating both Ttr and T„i at other pressures, we 
obtain the phase diagram of dense P-H2 in the T—P plane 
as shown in Fig. ^(a). The new solid phase is labeled as 
phase r. Since 0{R) vanishes, its structure is orientation- 
ally disordered. The I'/I and liquid/F phase boundaries 
are indicated by solid triangles and open circles with er- 
ror bars, respectively. Large curvature of the phase line 
for P < 50GPa reflects the rapid change of the volume in 
that region. A Simon equation of F = 1.2 x lO^'^T^'^^s 
(the dotted curves in Fig. ^) fits our results rather well. 
A lack of experiment on T,„ at high pressures makes a 
precise comparison between theory and experiment diffi- 
cult Jaut we obtain a reasonably good agreement with the 
dataEJ for P < 4GPa, as csm be seen in Fig. 1(b). The 
data for n-H2 up to 15GPa,Eil are always higher than our 
results, reflecting the fact that the presence of 0-H2 makes 
Tm higher In Fig. |^(a) we show our result for dense o- 
D2 exhibiting a similar phase diagram. The difference in 



Ttr between these two systems becomes smaller with the 
increase of P, but P-H2 gives definitely a higher Tfj. than 
0-D2. As for T,„, we find no meaningful difference. 



1000- 



800 



600 



400 



200 



(a) 0-D2 

liquid 



r {Cmca) 




I (hep) 



II 





50 



III 



(b) Classical 



liquid 



r {Cmca) 

with r|=0.61 J 




with r|=1 .0 



I (hep) 



150 



100 150 20050 100 

P (GPa) P (GPa) 

FIG. 5. (a) Similar phase diagram for 0-D2 with symbols 
indicating just the same meanings as those in Fig. ^a) . Two 
sets of experiments have detected three distinct phases I, II, 
and III for P below 200GPa at low temperatures, as indi- 
cated by two dashed and dotted-dashed curves with different 
thickness, taken from Refs. ^ and ^. In (b), solid triangles are 
the corresponding high-temperature phase transition points 
for the classical reference system (M — 1), while open squares 
show the transition points by the pair potential with rj = 0.61, 
defined in Eq. (||). (See the text in Sec. 0.) For comparison's 
sake, the double-dotted-dashed curve shows the data for solid 
H2, the same as in Fig. U (a). 



IV. QUANTUM HARDENING EFFECT 

For understanding the mechanism to bring about the 
phase F, we have performed another simulation with 
M = 1 (the "classical reference" system in which the 
molecules are treated as classical particles) and found an 
analogous phase diagram as shown in Fig. ^(b). In this 
system, only the temperature-induced random motion of 
molecules can play a role in transforming a compact hep 
structure into a less compact Cmca one and finally into 
a liquid. Thus we can identify the thermal fiuctuations 
as the basic driving force of the present phase transition. 

Examining the changes of various physical quantities at 
T around Ttr-, we found that the most important change 
occurs in the potential energy. The energy increases with 
T in both the classical (Fig. |(a)) and the H2 (Fig. |(b)) 
systems due primarily to the fiuctuation of molecules to 
shorten the nearest-neighbor distance Rq ^ 1. 9 A, be- 
cause V^air is very repulsive at R near i?o and it changes 
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very rapidly with R. (See Fig. |l|.) We should note that 
the increasing rate is larger in the compact hep phase 
than that in the less compact Cmca phase. In this re- 
gard, the phase transition is induced by the potential- 
energy gain as indicated in Fig. ^ No such gain is seen 
for the translational kinetic energy as shown in the inset 
ofFig.|(b). 
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FIG. 6. Temperature dependence of the potential energy 
(l^nucicar) per molecule for both (a) the classical reference sys- 
tem and (b) P-H2 near respective Ttr's at lOOGPa. Data in 
hep and Cmca structures are shown by solid and open sym- 
bols, respectively. The inset in (b) plots the change in the 
translational kinetic energy per molecule. 

Eventually the reason why Tf,. in the classical system 
is much lower than that in solid P-H2 must be ascribed 
to the quantum zero-point motion of protons, the only 
ingredient that is not included in the classical system. 
Such a motion is, in a sense, regarded as a kind of vac- 
uum polarization the effect of which, in the context of 
high-energy physics, can be renormalized into fundamen- 
tal physical quantities such as mass and interaction con- 
stant. We envisage that this must also be the case here. 
In order to substantiate this, we consider a toy model, 
namely, an interacting two-particle system in a harmonic 
potential in one dimension, described by the Hamiltonian 
as 



H = 



2m 



l^ + ^ + n.r-..), (4) 



from which we can derive the Hamiltonian for the relative 
motion as 



Hr = ^+ ^lUJ^X^ +Vix), 



(5) 



where /j, = m/2,uj = and x = xi — X2. By 

writing ^pi^) — '4'oi^)^{x) Stiid E = lo/2 + e with 
ipoix) — ^/JIuj/tt ex.p{—^ujx'^ /2), we can cast Eq. (|^) into 



P_ 
2^1 



+ yix) fix) = ef{x), 



(6) 



where the effective interaction V{x) is related to the bare 
one through 



V{x) = V{x) - 11-^^^ = V{x) 
Wo f 



(7) 



We see that the effect of the zero-point oscillation is in- 
cluded effectively in V{x), implying that the physics can 
be captured in terms of y(x).E3 In this sense, Vpair can 
include the effect of intramolecular vibrons. 

In general the closed-packed structure is known to be- 



come more stable with a harder Vpair B In the above toy 



model, we readily know that, if V{x) is repulsive, xip ftp 
is positive, leading V{x) to be more repulsive than V{x). 
This indicates that the zero-point motion of protons can 
make Vpair effectively more repulsive than that in the 
classical system in which V{x) is reduced to V{x) owing 
as TO —> 00. 
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R (A) 

FIG. 7. Schematic view to understand the Fluctuation ef- 
fect on the potential energy. Inset: (a) average potential en- 
ergy as a function of fluctuation 5R with R = 1.9 A. (b) 
average potential energy as a function of separation R with 
SR = 0.1 A. 

This hardening effect due to the molecular fluctuations 
can be understood more visually through Fig. |^. At the 
pressure range considered in this paper, the repulsive part 
of the interaction plays a central role in determining the 
potential energy. The fluctuation shortens one distance 
between the neighboring molecules, but at the same time 
it lengthens another one. Thus we have to average these 
two opposite effects. Because the potential changes very 
rapidly with R the intermolecular distance, AVi the in- 
crease of the potential due to the shortening is larger than 
AV2 the decrease due to the lengthening. Therefore the 
average potential energy, [V{R - SR) + V{R + 6R)]/2, 
becomes larger than V{R) the potential energy without 
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the fluctuation. As shown in the inset (a) of Fig. 0, the 
larger the fluctuation, the more repulsive the potential 
energy. This explains the behavior of the potential en- 
ergy in Fig. H; the solid H2 has a larger potential energy 
than the classical reference system at the same pressure 
and temperature, because it provides larger fluctuations 
due to the combination between the quantum and ther- 
mal effects. This also manifests itself in the much larger 
potential-energy gain in Fig. |(b) than that in Fig. |(a). 
Thus we can conclude that the molecules in the quantum 
system feel an effectively more repulsive potential than 
those in the classical system. This "quantum hardening 
effect" of Vpair explains why hep becomes more stable in 

With the increase of pressure, the deviation of the I/I' 
phase lines between solid H2 and the classical system is 
getting larger. In order to understand this result, we plot 
[V{R - SR) + V{R + SR)]/2 as a function of R with a 
given 6R in the inset (b) of Fig. 0. The average potential 
energy increases with the decrease of i?, a behavior of R 
seen with the increase of pressure. Thus we conclude that 
the quantum hardening effect becomes stronger with in- 
creasing pressure and this explains the above result well. 

In order to study this quantum hardening effect more 
quantitatively, we have performed another simulation in 
the classical system with the pair potential, defined as 

^pair = V^SG + VVSR + Vam, (8) 

where the parameter 77 is introduced to measure the soft- 
ness of Vpair; ?y = for the hardest SG potential and rj = 1 
for the softest SGH one. In Fig. |(b), the open squares 
with error bars show the I/F phase boundary obtained 
using V^air with 77 — 0.61. Notice that this more repulsive 
potential than the SGH one makes the I/F phase bound- 
ary of the reference classical system almost the same as 
the one for the quantum solid H2 with the SGH potential 
(the double-dotted-dashed curve) at P = lOOGPa. We 
need a softer (harder) V^air to better fit the data at lower 
(higher) pressure, but this behavior is consistent with 
the above-mentioned pressure dependence of the quan- 
tum hardening effect. 

V. SUMMARY 

Utilizing an intermolecular interaction potential Vpair 
proved to be reliable at high pressures in the PIMC sim- 
ulations, we have found the structural phase transition 
from orientationally disordered hep to Cmca induced by 
thermal fluctuations in both dense solid H2 and D2 be- 
fore melting. The result shows that a new phase exists at 
high temperatures in the solid molecular systems. The 
potential-energy gain due to the enhancement of fluctu- 
ations with increasing temperature is the main source to 
bring about the structural phase transition. The strong 
quantum zero-point motion plays an important role in 
determining Ttr the transition temperature through the 



quantum hardening effect on the potential energy. The 
melting temperature is also estimated for full under- 
standing the behavior of the solids and for the evaluating 
the temperature range of the new phase. 

Through additional simulation studies, we have also 
examined to see how our conclusion about the existence 
of the new phase is robust against the choice of Ipair- 
We have found that both Ttr and Tm are rather sensi- 
tive to the change of Vpair- An example of Ttr is given 
in Fig. 11(b). However, there is always a melting transi- 
tion after the I/F transition with increasing temperature, 
suggesting the robustness of the existence of the Cmca 
phase. 
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